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1.  Summary 

We  have  performed  under  the  subject  contact  a  theoretical  study  i)  to 
establish  a  mathematical  means  to  compute  the  electromagnetic  response  of  an 
earth  having  continuous  and  arbitrary  conductivity  variations,  and  ii)  to 
formulate  a  transfer  function  that  can  convert  a  wideband  frequency-domain 
response  to  a  continuous  conductivity-depth  function. 

The  study  is  based  upon  our  past  experiences  in  both  theoretical  and 
experimental  aspects  of  frequency-domain  electromagnetic  method.  We  present  a 
new  formal  integral  representation  for  Maxwell's  equations  for  a  medium  having 
an  arbitrary  conductivity  function.  It  is  shown  that  the  representation  is 
suitable  for  forward  problems  and,  through  a  proper  integral  transform,  the 
inverse  may  possibly  exist.  An  attempt  is  being  made  to  derive  a  transfer 
function  that  can  convert  an  electromagnetic  spectral  response  to  a 
conductivity-depth  profile. 


2.  Introduction  and  Background 

The  frequency-domain  electromagnetic  (EM)  method  employing  a  single 
frequency  or  several  frequencies  up  to  the  RF  band  has  long  been  in  use  in  the 
geophysical  industry  for  investigating  shallow  geological  structures  in  terms 
of  anomalous  electrical  conductivity.  The  method,  often  referred  to  as  the 
induction  EM  method,  involves  the  propagation  of  low-frequency  EM  waves  over 
the  earth. 

Ideally,  an  exploration  method  not  only  locates  a  sought  target  but  also 
is  able  to  describe  it  as  a  function  of  depth.  One  obvious  approach  which  can 
provide  both  penetration  and  vertical  resolution  is  a  method  employing  a 
wide-band,  multifrequency,  or  sweep-frequency  source. 

The  groundwork  for  the  EM  method  as  an  exploration  tool  for  conductive 
mineral  deposits  has  been  well  established  as  documented  by  Wait  (1962), 

Keller  and  Frischknect  (1966),  Ward  (1967),  Patra  and  Mallick  (1980),  Kaufman 
and  Keller  (1983)  and  others.  Although  the  method  is  based  on  the 
well-founded  classical  EM  theory,  the  applications  of  the  theory  to  realistic 
geological  conditions  have  been  limited  to  only  simple  cases. 

During  the  past  years,  the  principal  investigator  has  been  involved  in 
studying  the  wideband  induction  EM  response  through  both  theory  and 
experimentation.  Using  an  integral  equation  representation  for  solving 
Maxwell's  equation  having  inhomogeneous  boundary  conditions  (Won  and  Kuo, 

1975a  and  1975b),  we  computed  two-dimensional  spectral  profiles  of  a  circular 
cylinder  buried  in  a  conductive  half-space  under  a  line  source  excitation. 

The  secondary  EM  responses  (without  the  primary  field)  are  computed  for  a 
frequency  range  from  10  Hz  to  40  KHz  along  the  surface.  Ensembles  of  such 
spectral  profiles  in  a  distance-frequency  space  are  shown  in  Figure  1  for 
three  depths  of  the  circular  cylinder  (Won,  1980). 


In  order  to  test  the  theoretical  results,  we  experimented  with  a  scaled 
laboratory  sweep-frequency  EM  system.  The  conductive  earth  was  simulated  with 
saline  water  while  the  target  was  made  of  graphite.  Figure  2  shows 
experimentally  obtained  spectral  responses  for  a  vertical  graphite  slab 
submerged  at  three  different  depths  (Won,  1980).  Admittedly,  the  two  models 
are  dissimilar:  a  two-dimensional  model  of  a  circular  cylinder  with  a  line 
source  cannot  be  compared  rigorously  with  a  vertical  dike  in  the  presence  of  a 
dipole  source. 

The  spectral  characteristics  for  both  cases  have  been  qualitatively  known 
and  are  physically  understandable.  Rather,  the  main  advantage  of  such 
wide-band  spectral  data  is  the  method  of  cross-sectionally  displaying  the 
entire  data  to  provide  some  intuitive  interpretations :  An  anomaly  caused  by  a 
large  and  deep  structure  will  have  a  broad  feature  in  space  and  will  appear 
toward  the  lower  end  of  the  spectrum,  while  an  anomaly  caused  by  a  small  and 
shallow  structure  will  have  a  narrow  feature  in  space  and  will  appear  toward 
the  higher  end  of  the  spectrum. 

Encouraged  by  the  results  from  the  theoretical  study  and  laboratory 
equipment,  we  developed  and  field-tested  a  first  prototype  sweep-frequency  EM 
system  mounted  on  a  truck  (Won,  1982  and  1983).  A  transmitter  and  receiver 
pair  were  used  in  a  horizontal  coplanar  configuration  employing 
logarithmically  sweeping  harmonic  waves  from  500  Hz  to  80  kHz.  The  sweep 
duration  was  about  2  seconds  at  a  maximum  power  of  about  2  kW.  Secondary 
field  amplitude  and  phase  spectra  were  measured,  digitized,  and  stacked 
several  times  to  enhance  the  signal-to-noise  ratio.  Data  were  collected  while 
the  truck  was  in  transit  (at  about  5-10  MPH).  The  entire  operation  was 

i 

I 

automated  by  an  on-board  computer.  i 
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Figure  3  shows  one  of  the  spectral  profiles  obtained  over  a  metamorphic 
terrain  overlain  by  sediments  of  variable  thicknesses  with  sporadic  outcrops 
of  graphite  schist.  The  known  cultural  objects  as  well  as  available  surface 
geological  data  are  shown  on  an  approximate  geological  section.  The  anomalies 
due  to  cultural  objects  such  as  bridges  and  culverts  are  easily  recognizeable 
because  their  high-frequency  anomalies  are  big  enough  to  wipe  out  the  entire 
spectral  section.  In  contrast,  geological  anomalies  are  laterally  and 
vertically  extended  showing  a  spatially  distributed  anomaly  pattern.  This 
particular  profile  was  made  at  about  6  MPH  by  stacking  five  sweeps  per 
station,  each  sweep  lasting  2  seconds  and  having  70  discrete  logarithmic 
frequency  steps. 

A  spectral  profile  displayed  in  a  frequency-distance  space  is  somewhat 
analogous  to  a  reflection  seismic  profile  displayed  in  a  time-distance  space. 
Although  the  seismic  diffraction  theory  is  far  from  complete  and  more 
complicated  than  the  EM  diffraction  theory,  seismic  reflection  data  are  often 
self-evident  in  terms  of  the  relative  geometry  of  the  reflecting  boundaries. 

It  is  not  to  say  that  the  conversion  of  a  frequency-distance  space  EM  data 
into  a  conductivity-depth  section  would  be  as  definite  as  the  conversion  of  a 
seismic  time-distance  section  into  a  depth-distance  section;  the  spectral  EM 
theory  must  be  developed  much  further  to  reach  that  stage. 

In  an  analogy  of  a  seismic  (time-distance)  section  which,  given  the 
knowledge  of  the  velocity  distribution  everywhere,  can  be  converted  to  a 
depth-section,  we  may  surmise  that  an  EM  spectral  (frequency-distance) 
section,  given  the  knowledge  of  the  conductivity  distribution  everywhere,  may 
be  converted  to  a  similar  depth  section.  This  report  focuses  on  finding  the 
theoretical  linkage  between  the  conductivity  distribution  and  the  geological 
cross-section  through  the  understanding  of  wideband  EM  induction  responses  of 


the  earth 


We  have  also  applied  successfully  an  airborne  frequency-domain  EM  method 
to  shallow  ocean  bathymetry  (Won  and  Smit,  1985a)  as  well  as  measuring 
conductivities  of  water  and  bottom  sediments  (Won  and  Smits,  1985b).  Figure  4 
shows  an  exemplary  bathymetric  profile  (solid  line)  compared  to  ground  truth 
data  (solid  circles)  obtained  in  Cape  Cod  Bay.  We  used  two  frequencies  at  385 
Hz  and  7200  Hz  both  in  horizontal  coplanar  configurations  at  a  height  of  about 
50  m  above  the  ocean.  The  sensor  was  towed  by  a  helicopter  at  a  speed  of 
about  90  knots.  By  analyzing  the  dual  frequency  data  using  a  new  inversion 
algorithm,  we  were  able  to  determine  simultaneously,  water  depth,  water 
conductivity ,  and  sediment  conductivity.  The  method  has  a  potential 
application  to  determining  the  ice-thickness  profile  in  the  polar  region  (Won 
and  Smits,  1950).  A  separate  comprehensive  report  on  this  subject  is  appended 


to  this  Final  Report 


3.  Theoretical  Approach 


One  of  the  major  shortcomings  in  the  existing  theories  is  that  the  earth  I 

model  employed  for  their  forward  calculation  is  a  layered  earth  with  a  finite 
number  of  layers,  each  having  a  distinct  electrical  conductivity.  These  can 
rarely  be  more  than  three  to  five  layers  in  order  to  avoid  numerical 
instability  during  inversion.  Such  a  layering  may  be  geologically  unrealistic 
in  most  cases  (since  there  is  not  a  sufficient  number  of  layers  to  account  for 
gradual  changes)  and  is  believed  to  be  the  main  source  of  instability  in 
inverting  real  data.  Therefore,  the  following  equations  are  designed  for  the 
forward  calculation  of  a  laterally  homogeneous  earth  with  a  continuous 
conductivity  variation  in  depth. 

The  fundamental  equations  for  the  magnetic  field  generated  by  a  vertical 
oscillating  dipole  located  at  the  surface  of  a  horizontally  stratified  earth 
have  been  given  by  Wait  (1955),  Kozulin  (1963),  and  Frischknecht  (1967). 

The  following  is  the  expression  for  the  EM  vector  potential  F  generated  in  a 
point  at  or  above  the  surface  of  the  earth  under  conditions  considered  by 
Kozulin: 

00 

F  - 1 . -  +J  R(A,d  ,o  ,«)  J0(Ar)dA}  .  (1) 

Z  o  4  U 

/r2+(h-z) 

The  components  of  the  vector  potential  in  the  horizontal  plane  are  zero,  and 
the  radiation  term  in  the  Maxwell  equation  has  been  neglected. 

The  transmitting  dipole  is  considered  to  be  located  at  (0,-h)  of  the 
cylindrical  coordinate  system  with  the  z-axis  directed  vertically  downward. 

The  point  at  which  the  vector  potential  is  computed  is  located  at  (r,-z). 

The  mathematical  notations  are  as  follows: 

d[  thicknesses  of  the  i-th  subsurface  layer, 

conductivities  of  the  i-th  subsurface  layer. 
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^  V  V 


h  transmitter  height  (negative  upward) 
z  receiver  height  (negative  upward) 

X  variable  of  integration, 

M  dipole  moment  of  the  transmitter, 

u)  angular  frequency  of  the  transmitted  signal, 
u  magnetic  permeability  for  free  space,  and 

Jq,  zeroth  order  Bessel  function. 

The  kernel  function  R( X)  can  be  computed  from  the  subsurface  layer 
parameters  and  the  frequency  using  a  recurrence  relation.  Setting  R( X)  = 
Rq  n(X)  in  which  the  first  suffix  represents  consideration  of  the  field  in 
space  above  the  ground  surface,  and  the  second  suffix  n  is  the  number  of 
subsurface  layers,  the  recurrence  relation  is  given  by  [using  the 
notation  of  Koefoed  et  al.  (1972)]: 


R,,  (X) 

(i-1) ,n 


(1-D.1+  e-Mf-z) 


V 

1  +  V 


„  n  tRf  (X)e"2dtVi 
(i-l),i  i,n 


(2) 


and 


^n,n(X)  0 

where 

Vi  =  (X2  +  ki2) 1/2, 

2 

k  =  iunio(z) ,  and 

i 

Vi,k  =  (Vi  -  vk)/(Vi  +  vk). 

The  expressions  for  the  components  of  the  magnetic  field  strength  can  be 
obtained  from  (1).  The  relation  between  the  magnetic  field  strength  and  the 


vector  potential  can  be  written  as 


and  consideration  of  the  field  in  the  space  above  the  ground  surface  implies 
2  2 

that  k  =  ky  =  impose  =  0»  reducing  (3)  to 


H  =  V  V*F.  (4) 
Recalling  that  the  components  of  F  in  the  horizontal  plane  are  zero,  we  may 
rewrite  (4)  as 


H 


(f-  F  ) 

3z  z 


or 


H  = 
z 


and  H 

r 


3* 

3z  3r 


F 

z 


(5) 


Applying  the  differentiations  to  (1)  and  then  setting  z  =  0,  we  have,  for 
the  components  of  the  magnetic  field  strength  at  a  point  at  the  surface  of  the 
earth,  the  expressions 

00  2 

H  =  —  — — r  +  j  1  R(  A,d,  ,  o.  ,u)Jo(  Xr)dX,  (6) 

Z  4lTrJ  4ff  J0  1  1 


where  Jq  and  are  Bessel  functions  of  the  first  kind.  The  first  term  in  (6) 
represents  the  primary  field  in  the  absence  of  the  earth  while  the  second  term 
denotes  the  secondary  field.  Only  the  secondary  field  exists  for  the  radial 
component . 

The  recurrence  formula  (2)  is  widely  used  for  computing  the  F.M  responses 
of  a  multilayered  earth.  The  infinite  integrals  shown  by  (6)  and  (7)  are 
often  computed  numerically  using  a  Hankel  transform  algorithm  or  digital 
filter  method  (Koefoed  et  al,  1973;  Anderson,  1979  and  1982;  Chave,  1983). 

The  main  shortcoming  is  that  these  methods  are  strictly  limited  to  discretely 


lavered  earth  models. 


4.  Mathematical  Formulation  for  an  Earth  having  a  Continuous  Conductivity 


Variation 

We  now  consider  a  transformation  of  the  kernel  function  R( X)  in  an 
integral  form  which  can  handle  a  continuously-varying  conductivity  function 
o(z).  Our  present  goal  is  to  express  the  magnetic  field  as  an  integral 
equation  containing  a  continuous  function  a(z).  Such  expressions  are  highly 
desirable  for  interpreting  our  swept-f requency  data. 

In  order  to  change  R( A)  to  a  differential  form,  we  first  let  the 
thickness  of  each  layer  d^  be  Az,  a  vanishingly  small  distance.  Under  this 
condition,  the  terms  appearing  in  (2)  can  be  shown  as  follows: 


and 


V(i-1) ,i 


| - i -  /\^+k^(z)  }  •  Az 

2A“z+k2(z) 


-1 


4{\  +k2(z)} 


Ak 2(z) 
Az 


Az 


(8) 


d^V^  =  Az  J X2+k  z( z ) 


(9) 


Here  we  assume  that  k(z)  is  continuous  and  its  first  derivative  exists.  This 
is  found  to  be  the  main  condition  for  departing  from  the  traditional  layered 
earth  model.  It  can  be  seen  that  while  the  two  numerator  terms  In  (2)  are  in 
the  first  order  of  Az,  their  product  appearing  in  the  denominator  is  in  the 
second  order  of  Az.  For  the  limiting  case,  therefore,  we  may  write  equation 
(2)  as 

R(  i—  1 )  ’n  (X)  =  V(  i  —  1 )  ,i  +  Ri,n(X)e  2dlVi  * 

(10) 

Expanding  this  and  using  its  own  recurrence  relationship,  we  find 


R(A)  =  R0>n(A)  -  V0>1  +  V1>2e-2dIVl  +  V^e  2(diVt+d2V2> 

.  -2(d  iV.+doVo+d-jV.)  ,  -2(d,V,+...+d_V_)  , 


As  we  let  A z  +  0,  we  may  express  (11)  in  an  integral  form: 


R(X)  -  ~  J  j 

0 


1  iX 

V  dz 


z 

exp  {-2]  V(S)d?}dz 
0 


(12) 


where  V(z)  =  /Az+kz(z) .  The  exponential  term  denotes  physically  a  two-way 
attenuation  in  depth  while  the  remaining  terms  can  be  interpreted  as  an 
inductive  coupling  within  the  earth  medium.  Since  the  exponential  terra  is 
dominant,  we  would  generally  expect  a  stable  and  convergent  integration. 

g 

The  secondary  vertical  magnetic  field  Hz(u>,r),  the  second  term  in  (6), 
can  now  be  written  as 

QO  ^  OO  2 

«®U,r)  =  ‘  V  Ju(Xr>  J0  vfiT  exp(-2JoV(OdqdzdA  •  (13) 


Changing  the  order  of  integration  between  z  and  A,  and  using  the  relationship 


1  iX 

V  dz 


1 

2{A2+k2(z)} 


dk  2(z) 
dz 


(14) 


we  may  rewrite  (13)  as 


Hs(io,r) 
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co  ~  oo 

Mi  dk2(z) i 
8nJ o  dz  J  o 


a2J0(  ^g) 

A2+k2(z) 


z 

exp {-2  Jo 


^A2+k2(  O  d?}dAdz 


(15) 


The  A-integral  can  be  evaluated  by  the  residue  theorem.  We  first  note  that 
the  integral  is  symmetrical  in  A;  thus,  its  range  can  be  expanded  to  negative 
infinity.  We  find  two  simple  poles  in  the  second  and  fourth  quadrants  at 
A  =  ±k.  Also,  the  integration  along  an  infinite  semi-circle  in  the 
upper  half  of  the  complex  plane  can  be  shown  to  be  zero.  The  integral, 


therefore,  can  be  represented  by  the  residue  at  the  second  quadrant  and,  thus 


Hz(u),r)  =  —  j  Qk2(z)J0(ikr)  ^---exp  {-2  /k2(z)-k2( O  d  C }  dz.  (16) 

Equation  (16)  represents  a  significant  improvement  over  (6)  in  the  sense 
that  it  is  the  first  expression  by  which  we  can  compute  the  EM  response  over  an 
earth  model  having  a  continuous  conductivity  function.  Once  <j(z)  and  its 
derivative  are  specified,  the  integral  can  be  evaluated  either  analytically  or 
numerically  depending  upon  the  complexity  of  the  conductivity  function.  Its 
convergence  is  facilitated  owing  to  the  exponential  attenuation  term. 

Our  ultimate  goal  is  to  invert  equation  (16)  so  that  o(z)  can  be 
determined  at  each  data  point  for  a  given  wideband  spectrum  H(w,r)  at  r=a,  a 
constant  transmitter-receiver  separation.  We  plan  to  continue  our 
investigation  on  possibilities  of  analytically  inverting  (16)  through  integral 
transform  techniques.  Otherwise,  we  may  be  forced  to  solve  (16)  numerically. 

A  popular  technique  for  numerically  solving  integral  equations  is  the  moment 
method  as  described  by  Harrington  (1968).  The  method,  in  essence,  is  an 
approximate  linearization  process  which  assumes  a  certain  functional  for  the 
unknown  in  a  whole  range  or  in  each  subsectional  range.  This  process  reduces 
the  integral  equation  into  a  set  of  matrix  equations  which  is  inverted 
numerically. 

Equations  (6),  (7)  and  (16)  can  be  also  viewed  as  convolution  integrals: 
it  has  been  well  known  that  the  Hankel  transform  integrals  can  be  transformed 
to  a  linear  convolution-type  integral  through  proper  substitutions  of 
variables  (Ghosh,  1971;  Koefoed  et  al.,  1972;  Anderson,  1979  and  1982).  If 
the  continuous  convolution  is  discretized,  the  resultant  discrete  linear 
convolution  becomes  simply  a  linear  digital  filter  of  a  finite  length. 

We  may  formally  treat  Equation  (16)  is  a  similar  way:  the  filter  acts  on 


a  kernel  function  representing  the  conductivity  distribution  0(2)  and  its 
vertical  derivative.  The  filtered  output  corresponds  to  the  spectral  response 
of  the  magnetic  field  H(ui).  Therefore,  if  we  can  find  a  inverse  filter  (or  a 
transfer  function),  we  will  be  able  to  convert  the  spectral  response  H(u>)  to 
the  conductivity  function  o(z).  Finding  an  inverse  of  a  given  linear  digital 
filter  is  a  highly  advanced  subject  of  deconvolution,  and  a  variety  of 
approaches  are  available  particularly  for  seismic  signal  processing  (e.g., 
Robinson,  1983). 

5.  Conclusion 

Despite  the  fact  that  the  layered  earth  model  has  been  studied  by  many 
researchers,  the  model  has  not  been  used  comprehensively  in  terms  of  its 
spectral  responses.  Other  than  for  a  low  frequency  range  mostly  using  natural 
sources  such  as  those  in  MT  survey,  the  EM  method  measuring  continuous 
spectral  response  rarely  has  been  considered  as  a  stratigraphic  mapping  tool. 
While  the  EM  method  is  used  mainly  for  locating  isolated  conductive  orebodies, 
there  is  a  strong  theoretical  possibility  that  the  method  can  be  expaned  for 
wide  use  utilizing  the  spectral  profiles  displayed  in  a  frequency-distance 
space. 

The  progress  of  the  subject  research  is  still  elementary  considering  the 
complex  geological  nature.  Yet,  the  task  is  believed  to  be  a  timely  and 
essential  step  in  enhancing  the  capability  and  further  development  of  a 
wideband  (ground  or  airborne)  electromagnetic  exploration  method  applicable  to 
geological,  geophysical,  and  geotechnical  mapping  of  complex  three-dimensional 


earth  structures 
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the  center  of  the  profile.  Any  single-frequency  profile  may  be  obtained  by 
selecting  a  specific  frequency  on  the  vertical  scale.  The  contour  Interval  Is 
5  dB  with  an  arbitrary  reference. 


Figure  4.  AEM  bathymetric  profile  for  Line  5021  along  with  radar 
profile.  Solid  line  represents  AEM  bathymetry;  while  solid  clrcL 
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ABSTRACT 


Airborne  Electromagnetic  Bathymetry  (AEMB)  is  a 
recently  developed  technique  for  determining  the  depth  of 
shallow  sea  water,  using  an  airborne  electromagnetic  method. 
The  transmitting  and  receiving  loops  are  encased  in  a  rigid 
tube  and  towed  behind  an  aircraft  at  an  altitude  of  about 
140  feet.  The  relatively  simple  geometry  of  the  shallow  sea 
water  environment  provides  us  some  practical  advantages  in 
survey  design  and  data  interpretation,  while  the  small  skin 
depth  associated  with  highly  conductive  sea  water  requires 
low  frequency  which  makes  the  AEMB  problem  a  ma;or  chal¬ 
lenge  . 

Application  of  gee prys i cal  inverse  theory  'nonlinear 
and  linear,  combined  witr.  two-layer  model  simulating  sea 
water  and  bottom  s-d  :.:~r.*  provides  us  valuable  information 
toward  survey  d~s.gr  an:  da*a  interpretation.  Uncertainty 
in  the  sensor  alti* is  resolved  tnrough  the  inversion 
process  using  high  frequency  data  for  which  the  ocean  may  be 
considered  as  an  infinite  naif-space.  The  effect  of  varia¬ 
tion  of  the  coil  axis  witr.  respect  to  the  sea  surface  can  be 
resolved  if  the  inclination  is  measured  and  incorporated 
into  the  inversion  process.  An  optimum  frequency  range  for 
a  given  bathymetric  range  may  be  determined  in  advance  by 
analyzing  the  information  density  matrix.  For  a  reliable 
bathymetric  determination,  at  least  one  low  frequency  whose 
skin  depth  is  close  to  water  depth  need  be  included.  High 


frequencies  whose  skin  depths  are  less  than  1/4  -  1/3  of  the 
water  depth  do  not  contribute  significantly  to  the  bathyme¬ 
tric  solution,  however  the  high  frequency  data  can  be  used 
to  determine  the  conductivity  of  sea  water  and  to  correct 
for  the  error  in  the  altimeter  height  of  the  aircraft. 

Since  the  conductivity  of  the  sediment  is  proportional 
to  the  porosity  of  sediments  following  Archie's  Law,  it  is 
expected  to  be  decaying  with  depth.  A  theoretical  electro¬ 
magnetic  response  over  exponentially  decaying  conductivity 
material  underlain  by  sea  water  has  been  derived  and  com¬ 
pared  with  that  of  layered  earth. 


INTRODUCTION 


Recently,  there  has  been  an  increasing  interest  in  pos¬ 
sible  airborne  determination  of  bathymetric  depth  over  shal¬ 
low  water  (lakes  and  coastal  bays).  Becker  et  al .  (1984) 
have  considered  the  Airborne  Electromagnetic  Bathymetry 
problem  by  time-domain  electromagnetic  (TDEM)  method.  In 
this  paper,  we  will  present  a  fundamental  theory  for  the 
Airborne  Bathymetry  by  electromagnetic  frequency  sounding 
method . 

Unlike  cases  of  buried  conductive  layers  under  resis¬ 
tive  overburden,  the  airborne  electromagnetic  bathymetry 
(AEMB)  is  to  determine  the  depth  of  highly  conductive  sea 
water  over  relatively  resistive  bottom  sediments  (Figure  1). 
It  is  understandably  a  difficult  problem  to  determine  the 
bottom  geometry  of  the  sea  floor  by  electromagnetic  induc¬ 
tion  method;  because  of  the  high  conductivity  of  sea  water, 
EM  energy  attenuates  rapidly  with  depth  and,  thus,  the 
effect  of  the  bottom  geometry  on  the  returned  signal  is 
small.  The  resolution  is  dependent  upon  the  accuracy  of  the 
measured  data  and  their  frequency  content. 

In  the  airborne  survey,  the  transmitting  and  receiving 
loops  are  usually  encased  in  a  rigid  tube  or  a  truss  support 
'generally  called  a  ’bird’)  which  is  towed  behind  the  air¬ 
craft  on  a  guywire.  Height  is  measured  by  a  radar  altimeter 
usually  mounted  on  the  aircraft.  Horizontal  levelling  of 
the  bird  could  be  monitored  by  an  inclinometer.  Modern 


instruments  under  ideal  conditions  can  resolve  the  secondary 
EM  response  up  to  1  ppm  or  better  with  respect  to  the  pri¬ 
mary  field.  A  study  in  this  paper  shows  that  small  errors 
in  the  positioning  of  the  bird  can  cause  as  much  effect  as 
that  of  bottom  geometry  of  sediments,  indicating  that  cor¬ 


rection  of  the  data  must  be  made  before  interpretation. 
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Because  of  the  skin  effect  of  EM  wave,  the  deeper  sea 
water  becomes,  the  lower  frequency  we  need  to  properly 
resolve  the  deep  structure.  Technically,  it  is  more  diffi¬ 
cult  to  generate  a  lower  frequency  (below  100  Hz)  than  to 
generate  a  higher  frequency;  the  large  dipole  moment 
demanded  for  low  frequency  requires  an  extra  large  coil  and 
a  high  current  source  which,  in  turn,  increases  on-board 
logistics,  weight,  and  power  requirements.  In  addition,  the 
spatial  resolution  also  suffers  as  frequency  is  lowered. 
Therefore  it  is  ideal  that  optimal  range  of  frequency  be 
determined  for  a  given  range  of  water  depth.  The  result  of 
this  study  shows  that,  to  determine  basement  geometry  with 
data  measured  at  10  m  altitude  and  10  m  separation,  we  need 
at  least  a  frequency  whose  skin  depth  is  the  same  as  the 
depth  of  investigation,  and  we  do  not  need  those  frequencies 
whose  skin  depths  are  less  than  1/3  -  1/4  of  depth  of  inves¬ 
tigation. 

Once  the  data  are  corrected  properly,  the  parameter 
most  strongly  influencing  the  data  is  the  lateral  variation 
of  conductivity  of  sea  water.  The  next  most  significant 
parameter  is  the  thickness  of  the  sea  water  layer.  The  con¬ 
ductivity  of  bottom  sediments  has  the  smallest  effect. 
According  to  Grant  and  West  (1965),  sea  water  conductivity 
varies  between  2.5  and  5*0  (mho/m),  and  normal  basement  con¬ 
ductivity  is  in  the  range  of  0.0  to  0.4  'mho/m).  The  con¬ 
ductivity  of  basement  water  or  transition-zone  sediments 
varies  between  0.0001  and  1.0  (mho/m).  Since  the  target  of 


an  AEMB  survey  may  be  considered  to  be  shallow,  increasingly 
deepening  water  (less  than  100  m  in  depth)  on  the  continen¬ 
tal  shelf,  we  will  simplify  the  model  as  follows: 

i)  the  stratification  of  conductivity  within  sea  water  will 
be  neglected, 

ii)  since  the  slope  of  sediments  is  gentle,  the  boundary  of 
water  and  sediments  will  be  assumed  to  be  horizontal  for 
each  set  of  data,  and 

iii)  the  conductivity  variation  of  sediments  with  depth  will 
be  neglected  until  later  section.  However,  to  account  for 
this  oversimplification,  upper  and  lower  bounds  will  be 
sought  for  the  solution. 

These  simplifications,  lead  us  to  deal  with  a  two-layer 
model  with  three  unknowns,  viz,  sea  water  conductivity  '  01  ) , 
water  depth  (d),  and  bottom  sediment  conductivity  Ok.  The 
interpretation  technique  employing  Marquardt's  (1963,  1970) 
technique  and  generalized  inverse  theory  '(Wiggins,  1972; 
Jackson,  1972;  Inman  et  al.,  1973)  can  be  used  for  this  two- 
layer  model.  Applied  to  synthetic  data,  the  technique  pro¬ 
duces  an  accurate  determination  of  the  sea  water  conductiv¬ 
ity.  The  determination  of  two  other  unknowns,  depth  of  sea 
water  and  conductivity  of  sediments,  is  dependent  upon  the 
frequency  range  used  and  errors  involved  in  the  data.  'Edge 
model'  is  defined  as  a  maximum  or  minimum  depth  of  water  for 
a  given  condition;  since  sea  water  conductivity  is  deter¬ 
mined  accurately,  the  lower  bound  is  found  by  expanding  the 
water  depth  until  the  rms  difference  between  the  data  and 
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the  predicted  response  becomes  a  tolerable  preset  value, 
assuming  an  insulating  ocean  bottom.  Similarly,  the  upper 
bound  is  found  by  assuming  0.4  mho/m  as  the  maximum  conduc¬ 
tivity  of  bottom  sediments. 

We  will  first  review  the  response  function  of  vertical 
and  horizontal  dipoles  over  a  layered  earth.  Next,  the 
efforts  will  be  concentrated  on  the  method  of  correcting  the 
errors  which  can  be  caused  by  mispositioning  of  dipoles  in 
airborne  measurements,  and  on  the  optimum  frequency  range 
design  of  airborne  EM  bathymetry  by  using  the  generalized 
inverse  theory.  An  application  to  synthetic  data  will  be 
illustrated.  Finally,  a  theoretical  electromagnetic 
response  over  exponentially  decaying  conductive  material 
underlain  by  sea  water  will  be  presented. 


RESPONSE  OF  A  WATER  LAYER  OVER  SEDIMENTS 


In  the  very  low  frequencies  (VLF)  where  the  frequency 
range  is  below  100  kHz,  the  effect  of  the  dielectric  con¬ 
stant  which  governs  the  displacement  current  within  the  tar¬ 
get  materials  is  small  (for  water  € /*.  =  80).  Magnetic  sus¬ 
ceptibilities  of  sea  water,  or  sediments  are  essentially 
zero,  the  same  as  that  of  air.  For  most  crystalline  rocks, 
it  is  small  unless  they  contain  considerable  amounts  of  mag¬ 
netite.  The  main  physical  property  governing  the  electro¬ 
magnetic  field  is  the  electrical  conductivity .  Therefore, 
the  response  functions  of  layered  earth  model  (Wait,  1951; 
Frischknect,  1967;  Ward,  1967;  Patra  and  Mallick,  1980; 
Kaufman  and  Keller,  1985)  are  valid  for  the  bathymetry  prob¬ 
lem  . 

Since  the  transmitting  and  receving  loops  are  encased 
in  a  rigid  tube  or  a  truss  support  in  the  AEMB  measurement, 
the  primary  field  is  dependent  upon  radial  distance  only. 
Using  similar  notation  as  Frischnecht  (1967),  the  magnetic 
fields  due  to  a  vertical  magnetic  dipole  over  a  layered 
earth  (Figure  1)  become 
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and  those  due  to  a  horizontal  magnetic  dipole, 
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where  Hz  is  the  z  component  of  a  primary  field  due  to  a  ver- 

p 

tical  magnetic  dipole,  Hr  is  the  r  component  of  a  primary 
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field  due  to  a  horizontal  magnetic  dipole,  and  Hz  and  Hr  are 
the  z  component  and  r  component  of  secondary  magnetic  field. 


respectively. 

M  is  the  dipole  moment,  r'  is  the 

distance 

from  the  transmitting 

dipole,  and  r  is  the  horizontal 

compo- 

nent  of  r ' . 
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where  A  is  a 

variable 

of  integration  and  R(A)  is 

the 

kernel 

function  related  to  the  geological  parameters.  The  h  and  z 
are  the  heights  of  transmitting  and  receiving  loops,  respec¬ 
tively.  Jo  and  J,  are  Bessel  functions  of  order  zero  and 
one,  respectively.  For  a  two-layer  model  the  kernel  func¬ 
tion  is  given  by  (Patra  and  Mallick,  1980) 
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where  nt  =  y  A*  +  jw/xff;  in  which  i  represents  each  layer,  ' 
yC7,  w  angular  frequency,  ju  magnetic  permeability  and  Ct  is 
the  conductivity  of  each  layer. 


Dividing  secondary  fields  by  their  primary  fields,  we 


obtain  mutual  impedance  ratios; 


Z,  = 

-  r*  t0  , 
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where  Z3  = 

-  (/'  /r  )  t*  • 
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Z,  represents  a  mutual  impedance  ratio  of  vertical-vertical 
dipoles.  Zz  represents  a  mutual  impedance  ratio  of  perpen¬ 
dicular  dipoles  (when  the  source  dipole  is  horizontal,  Z* 
must  be  multiplied  by  two) .  Z^  represents  a  mutual  impe¬ 
dance  ratio  of  horizontal-horizontal  dipoles  whose  axes  are 
in  the  same  direction,  whereas  Zj  represents  that  of  hori¬ 
zontal-horizontal  dipoles  which  are  in  the  common  plane. 
When  the  heights  of  transmitting  and  receiving  loops  are 
kept  the  same,  the  configurations  are  called  horizontal 
coplanar,  perpendicular,  vertical  coaxial  and  vertical 
coplanar  systems,  respectively  (Figure  2). 


Transmitter 

Receiver 

-0- 

O 

O 

$ 

Horizontal  coplanar  (Z,) 
Perpendicular  (Z2) 
Vertical  coplanar  ( Zj  ) 
Vertical  coaxial  ( Z«. ) 


Figure  2  Four  different  configurations  for  vertical  and 
horizontal  magnetic  dipoles. 


Several  Z,  response  curves  are  shown  in  Figures  3 
through  5  for  a  two-layer  model  of  which  the  first  layer  is 
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a  highly  conductive  layer  simulating  a  sea-water  layer 
underlain  by  a  uniform  half-space  of  resistive  bottom  sedi¬ 
ments.  The  O;  and  (J2  represent  the  conductivity  of  water  and 
sediment,  respectively.  The  d  represents  thickness  of  the 
water  body.  In  Figure  3,  quadrature( imaginary) , 
inphase ' real)  and  amplitude  of  the  Z{  response,  in  ppm  are 
illustrated  for  a  two-layer  model  (top)  of  which  CT,  =  2 
mho/m,  (Jv  =  0.2  mho/m  and  d  =  12  m.  The  response  curves  are 
computed  at  a  10  m  height  with  a  10  m  coil  separation  in  a 
frequency  range  of  1 0  Hz  and  100  kHz,  and  presented  in  log- 
log  scale.  Figure  4. a  shows  the  effect  of  changing  depth 
while  keeping  the  conductivities  and  other  parameters  the 
same  in  the  model.  Figure  4-b  shows  the  effect  of  changing 
CJ,  •  Figure  4-c  shows  the  effect  of  changing  Oi. .  Figure  5 
illustrates  the  effect  of  changing  r  and  h.  Figures  4. a, 
4.b,  and  4-c  show  clearly  that  for  the  fixed  r  and  h,  the 
conductivity  0)  of  first  layer  has  the  most  significant 
effect,  followed  by  depth,  and  O'*  in  a  descending  order  in 
signif icance .  One  interesting  feature  is  that  the  skin 
depth  of  water  at  a  frequency  where  the  imaginary  component 
becomes  a  maximum  coincides  roughly  with  the  depth  of  water. 
In  Figure  5,  we  also  notice  that  inaccuracy  in  height  and 
separation  of  the  bird  can  cause  significant  effect  for  the 
geological  interpretation. 
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Figure  3  Mutual  impedance  ratio  of  a  two-layer  model 

'top),  in  ppm.  Amp,  Re  and  Im  represent  Amplitude,  Real 
and  Imaginary  components  of  horizontal  coplanar  system, 
respectively . 


Figure  4  Effects  of  parametric  changes  for  the  model 

shown  on  Figure  3  'O’,  =  2  mho/m,  Oi  =  0.2  mho/m,  d  =  12  m 
r  =  h  =  10  m),  while  keeping  the  other  parameters  the  same 

a.  Effect  of  changing  depth  (5,  10,  20,  40,  80  m) 

b.  Effect  of  changing  water  conductivity  '10,  5,  2,  1, 

0.5  mho/m) 

c.  Effect  of  changing  bottom  conductivity  '2,  1,  0.4, 

0.1,  0.001  mho/m) 
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CORRECTIONS  FOR  MISPOSITION  OP  DIPOLES 


In  the  airborne  measurement,  there  are  three  possible 
errors  in  positioning  the  transmitting  and  receiving  loops; 
i)  deviation  of  the  bird  from  a  horizontal  plane,  ii)  change 
of  the  distance  between  the  two  loops,  and  iii)  error  in  the 
altimeter  height  of  the  bird.  This  section  will  be  devoted 
to  investigating  the  effect  of  each  of  these  errors  in  the 
horizontal  coplanar  and  vertical  coaxial  systems.  We  will 
also  discuss  methods  of  correcting  for  these  errors. 

Effect  of  The  Bird  Inclination 

Let  us  consider  a  rigid  bar  along  which  horizontal 
coplanar  coils  are  mounted.  The  bar  is  inclined  at  an  angle 
9  with  respect  to  the  surface  as  shown  in  Figure  6. 


Figure  6 


Horizontal  coplanar  coils  tilted  by  an  angle  0. 
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As  long  as  the  tilt  angle  for  transmitting  and  receiving 
loops  remains  the  same,  the  primary  field  does  not  change. 
Verma  (1973)  have  considered  the  effect  of  ground  slope  in 
combination  with  transmitter  coil  misorientation ,  neglecting 
the  effect  on  the  primary  field.  The  secondary  fields  are 
not  affected  by  the  heights  of  trasmitting  and  receiving 
loops  individually,  as  long  as  the  center  of  two  receivers 
are  kept  the  same,  i.e.,  (h+z)/2  =  const.  However,  in  this 
case  the  horizontal  distance  between  the  two  receivers  is 
changed  f  r  — >  r*  cos9) .  Therefore,  the  error  can  arise 
from  the  angular  deviation  9  and  the  change  in  the  distance 
r  . 


The  magnetic  dipole  moment  of  the  transmitting  loop  can 
be  decomposed  by  M  =  M  cos9  z  +  M  sin9  r  .  This  decomposi¬ 
tion  results  in  the  same  effect  as  having  two  transmitting 
dipoles:  a  vertical  dipole  having  a  moment  of  M  cos9,  and  a 

horizontal  dipole  having  a  moment  of  M  sinB.  The  magnetic 

VS  VS  US  MS 

field  at  receiving  loop  becomes  H  — >  H  cos9,  H  — >  H  sin9, 
where  the  complementary  superscripts  V  and  H  are  to  denote 
vertical  and  horizontal  source  dipole  origin,  respectively. 
On  the  receiving  loop,  we  measure  the  field  as  a  vectorial 
sum  of  the  horizontal  and  vertical  fields  ;'HS  =  cos9  +  H^ 
sin9).  Therefore  the  expression  for  the  secondary  field 
measured  at  the  tilted  receiving  loop  becomes: 


s 

H 


'(cos©,  ,  sin9,  ( 

r 

Hz 

Vi  ' 

Hr 

’  c  o  s  6Z 

MS 

ns 

L  Hz 

Hr 

sinG, 

'15' 


where  Hz  :  z  component  of  secondary  magnetic  field  caused  by 


vertical  source  dipole, 

HS 

Hz  :  z  component  of  secondary  magnetic  field  caused  by 
horizontal  source  dipole, 

VS 

Hr  :  r  component  of  secondary  magnetic  field  caused  by 
vertical  source  dipole,  and 

HS 

Hr  :  r  component  of  secondary  magnetic  field  caused  by 
horizontal  source  dipole. 


0,  and  0i denote  the  inclination  angle  of  the  transmitting  and 
receiving  dipoles,  respectively.  The  first  bracket  on  the 
right-hand  side  of  equation  (15)  represents  the  effect  of 
tilting  the  source  dipole  by  an  angle  0, ,  while  the  third 
represents  the  effect  of  tilting  the  receiver  dipole  by  an 
angle  6*.  The  expression  (15)  can  be  verified  by  consider¬ 
ing  following  four  different  settings. 


i )  e,  =  o’, 

02  =  0 

4> 

i  i )  0,  =  0°, 

=  90* 

-0- 

i  i  i  )  0,  =  90*, 

02  =  o* 

-o- 

c> 

iv)  6,  =  90*, 

01=  90* 

~0~ 

From  equations  ' 2 ) 

through 

'  6 )  we 

have , 

H 


vs 

Hz 


Vi 

Hr 


S  HS 

H  =  Hz 


hS 

=  Hr 


HS 

Hz 


oHS 

dr 


vS 

-  Hr 


Hit 

HU.  r  z 


'16; 

17' 


Substituting  equations  (16)  and  '17)  into  15'  and  letting  0( 
=0=0  ,  we  obtain  the  expression  for  the  secondary  mae- 


netic  field 


HS  =  Hz  -  {--  -  tr)  sin1©  '18) 

47L 

Corresponding  mutual  coupling  ratio  Z”  is  obtained  by  divid¬ 
ing  equation  (18)  by  a  primary  field  - 

z ;•  (=  hs/hz)  =  z;  -  z ;  sinae  (19) 

where  Z\  -  -  r*  t0(r,cos9)  '20) 

Z'  =  -  r*  t2,(r'cos0)  /  cosO 

Similarly,  for  the  horizontal  dipole  (Figure  7)  the  secon¬ 
dary  field  measured  at  tilted  receiver  dipole  becomes 


s 

H 

[cos©,  ,  -sin©,  ] 

* 

HS 

Hr 

> 

HS 

Hz 

*  > 
COS0Z 

i/>  f 

w 

Hz5 

> 

-sinGz. 
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Figure  7  Vertical  coaxial  coils  tilted  by  an  angle  ©• 


and  corresponding  mutual  impedance  ratio  Z£  is 


z-  (=  HS/wPr)  = 

z;  -  '1  /2)  •  z;  sin1  0 

'22) 

where  Z^  =  ( 1 /2  )  • 

'K  -  z.'> 

'23) 

The  Z('  ,  Z3  and  Z'f  given  by  equation  '  20)  and  '23)  are  func¬ 
tions  of  r'cos©,  while  the  2,  ,  Z3  and  Zy  in  equations  (11)  - 
(14)  are  functions  of  r  (  =  r'). 

The  first  terms  on  the  right-hand  side  in  equation  '19) 
and  (22)  represent  the  effect  of  a  change  in  the  horizontal 
distance  between  two  loops  caused  by  tilting,  while  the  sec¬ 
ond  terms  are  the  effect  of  tilt  angle. 

Table  l.a  shows  the  real  part  of  Z,  ,  Z ,  -  Z^  sin20, 
and  Z"  of  tilted  horizontal  coplanar  system  for  0=5°.  The 
response  function  is  the  mutual  impedance  ratio  in  ppm,  for 
a  two-layer  model  where  r  =  10m,  z  =  h  =  10  m,  (T,  =  2  mho/m, 
d  =  12  m,  and  (7*.=  0.2  mho/m  (the  same  as  shown  in  Figure  3). 
In  the  table,  the  second  column  ( Z,  )  represents  the  response 
of  a  horizontal  coplanar  loop  at  corresponding  frequencies 
(given  in  the  first  column)  for  a  horizontal  bird.  The 
third  column  (Z/ )  shows  the  response  with  changed  distance 
(r'  =  r  cos©).  The  fourth  column  (-Z3sinz0)  shows  the 
effect  of  tilt  angle  (0  =  5°)>  and  the  fifth  column  '  Z " )  the 
total  expected  response  (sum  of  the  third  and  fourth  col¬ 
umns)  .  The  numbers  in  parentheses  in  the  third  and  fifth 
columns  represent  the  deviation  from  the  response  of  a  hori¬ 
zontal  coplanar  system  ( Z, ) - 


Because  the  tilt  angle  can  be  measured  by  an 
inclinometer,  the  horizontal  distance  (r'cos©)  can  be  calcu¬ 
lated.  Therefore,  if  we  use  r'cos©  as  a  distance  instead  of 
r  ,  the  effect  of  change  of  distance  is  not  problematic.  In 
this  model  ((T,  =  2  mho/m,  01=  0.2  mho/m,  d  =  12  m,  r  =  10  m, 
(h  +  Z)/2  =  10  m) ,  the  actual  change  in  distance  caused  by  a 

o  .  •  . 

5  tilt  angle  is  3*8  cm  1 0  m  -  1 0  m  x  cos5  )•  The  paren¬ 
thetical  numbers  in  the  third  column  of  Table  1 .a  show  the 
effect  of  a  3-8  cm  change  in  distance. 

The  second  column  in  Table  1 .b  lists  the  resultant 
field,  Z" ,  due  to  the  tilt.  This  term  is  considered  to  be 
measured  data.  The  third  column  in  the  table  shows  the 
effect  of  tilt  angle  for  a  half-space  (infinite  water 
depth) .  Comparison  of  this  column  with  the  fourth  column  of 
Table  l.a,  which  is  the  effect  of  a  tilt  angle  over  12  m 
deep  sea  water,  shows  that  except  for  small  deviations  at 
intermediate  frequency  ranges  (100  -  1000  Hz),  the  differ¬ 
ence  in  amplitude  is  trivial.  The  contribution  of  water 
depth  to  a  tilt  angle  effect  is  small.  Therefore,  the  cor¬ 
rection  for  the  effect  of  tilt  angle  can  be  made  by  sub¬ 
tracting  the  theoretically  calculated  effect  over  a  half¬ 
space.  The  last  column  in  Table  1 .b  was  obtained  by 
subtracting  the  third  column  from  the  second  column.  The 
numbers  in  the  parenthesis  represent  the  deviation  of  this 
column  from  the  Z\  response  in  the  third  column  of  Table 
l.a.  If  we  use  a  rough  estimate  of  depth  in  calculating 
'  -Z3’  sin1  9)  ,  we  will  have  more  accurate  Z'  response  and 


decrease  the  deviation.  Since  the  tilt  angle  8  can  be  meas¬ 
ured,  we  will  use  this  Z('  (r’cos0)  as  a  repsonse  function 
for  the  inversion.  If  there  are  other  factors  affecting  the 
change  of  separation  of  two  loops  other  than  the  tilt  angle 
(e.g.,  thermal  expansion  of  the  solid  bar  etc.),  those  also 
must  be  considered  in  calculating  Z('  response. 

Imaginary  component  of  horizontal  coplanar  system,  and 
real  and  imaginary  component  of  vertical  coaxial  system  also 
show  similar  relations  (Son,  1985).  If  the  tilt  angle  0  is 
small,  both  the  effect  of  tilt  angle  (-Z^ sin* 0)  and  the 
change  in  distance  r'cos©  become  small.  Depending  on  the 
amount  of  tilt  angle  (measued  by  an  inclinometer),  it  is 
necessary  to  determine  whether  the  correction  must  be  made 
or  can  be  neglected.  The  tolerable  range  of  0  is  highly 
dependent  upon  the  conductivity  of  the  water. 


Table  l.a  Effect  of  5  tilt  angle,  real  component  of 
horizontal  coplanar  system,  in  ppm. 


f  (Hz) 

z, 

2; 

-ZjSin*  6 

zv 

1 

1 

1  '0) 

0 

1  '0) 

10 

41 

41  '0) 

-0 

41  '0) 

100 

1978 

1979  (D 

-8 

1971  (-7) 

1000 

32132 

32172  (40) 

-132 

32039  '-93) 

10000 

86829 

87105  '216) 

-412 

86693  ' —1 36 ) 

100000 

113120 

113610  (490) 

-590 

113020  (-100) 

Z,  : 

no  tilting 

(  r  =  r  * ) 

r 1 cosQ) ; 

z  •  • 

change  in 

separation  (r  = 

values  in  pare 

nthesis  denote  deviation  from  Z,  due  to  change  in 
separation 

-Z3sin*  0  :  error  due  to  tilting  (0  =  5*  ) 

Z"  :  resultant  field  due  to  the  tilt;  values  in  paren¬ 
thesis  denote  total  deviation,  i.e.,  Z,  -  Z" 


Table  l.b  Correction  for  the  effect  of  5  tilt  angle, 
real  component  of  horizontal  coplanar  system,  in  ppm. 


f  (Hz) 

zr 

^-Z3' sin*  Q)l/z 

z"  +  'z; 

sinz 

a> 

> 

1 

i 

-0 

i 

(0) 

10 

41 

-0 

41 

(0) 

100 

1971 

-15 

1  986 

'!) 

1000 

32039 

-134 

32173 

(1) 

10000 

86693 

-412 

871  05 

'0) 

100000 

113020 

-590 

1 13610 

'0) 

’  —  Z j  sin1 

0 )yt :  error 

due  to  tilting 

'  0  =  5°  )  over 

half- 

space 

Z"  +  ' Z3’ sirf  0 )y •  :  correction  by  subtracting  half-space 

effect  from  Z"  . 

Values  in  the  parenthesis  represent  deviation  of  corrected 
term  from  Z! • 
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Effect  of  Error  in  Altitude 

In  this  section  we  will  look  at  the  effect  of  error  in 
altitude  of  the  bird  and  discuss  how  to  correct  for  the 
error.  Table  2. a  shows  the  effect  of  a  1  #  altitude  error 
(10  cm)  in  positioning  of  height  for  the  same  model  shown  on 
Figure  3 •  The  second  column  in  the  table  shows  the  real 
component  of  the  difference  between  Z,  ■  r  =  10  m)  and  Z''r  = 
9.9  m) .  This  table  shows  that  a  1  £  error  in  positioning 
height  is  more  problematic  than  a  5°  bird  inclination. 

Referring  back  to  Figure  4. a,  we  note  that  at  high  fre¬ 
quencies  the  response  converges  asymtotically  regardless  of 
the  depth  of  water.  Table  2.b  shows  how  the  response  func¬ 
tion  converges  as  a  function  of  depth  and  frequency.  The 
response  becomes  identical  up  to  1  ppm  for  water  depth 
beyond  20  m  at  10  kHz,  12  m  at  25  kHz,  and  8  m  at  50  kHz. 
This  saturation  in  the  high-frequency  response  will  occur  at 
a  shallow  depth  as  either  or  h  increases.  A  similar  phe¬ 
nomenon  occurs  in  the  response  attributable  to  the  sediment 
conductivity.  This  fact  renders  us  to  use  the  high  fre¬ 
quency  data  to  invert  the  height  of  the  bird  together  with 
the  water  conductivity.  For  example,  the  real  and  imaginary 
components  of  Z,  and  Z^  at  frequencies  25  and  50  kHz  '8 
data)  can  be  used  to  determine  the  height  of  the  bird  and 
conductivity  of  sea  water,  provided  that  the  sea  water  is 
deeper  than  12  m.  Table  3*a  illustrates  the  data  contami¬ 
nated  by  a  0.5  iS  random  error,  and  that  without  error.  The 
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data  with  0.5  %  error  were  subjected  to  the  nonlinear 
least-squares  (NIS)  inversion  to  determine  the  height  of  the 
bird  and  sea  water  conductivity,  the  theoretically  fitted 
values  are  listed  in  the  last  two  columns  of  the  Table.  The 
inverted  result  are  shown  in  Table  3*b.  For  initial  values 
of  h  =  8  m,  (X,  =  1.6  mho/m,  the  inversion  routine  finds 
almost  the  true  answer  h  =  10  m,  (T,  =  2  mho/m.  In  actual 
measurements,  the  particular  frequency  must  be  determined 
experimentally  for  a  given  bathymetric  range. 


Table  2. a  Effect  of  a  1  £  altitude  error  '10  cm)  at 


an  altitude 

of  10  m . 

Deviations 

from  correct 

response 

are  listed 
'  h  =  z  =  r 

in  ppm. 

=  1 0  m,  0( 

=  2  mho/m, 

01  =  0.2  mho/m, 

d  =  12 

f  (Hz) 

Re' Z,  ) 

Im( Z,  ) 

Re'Z4) 

Im '  Z^) 

1 

0 

-1 

0 

0 

10 

0 

-1  1 

0 

+2 

100 

-13 

-102 

3 

20 

1000 

-406 

-548 

-97 

87 

10000 

-1490 

-536 

1  96 

-9 

1 00000 

-1  960 

-1  83 

99 

-54 

Table  2.b  Responses  at  high  frequency,  in  ppm. 

The  underlines  indicate  depth  at  which  the  saturation 


occurs 
( h  =  z 

tor  each  frequency. 

=  r  =  10  m,  O’,  =  2  mho/m,  01  = 

0.2  mho/ 

m,  d  = 

12  m). 

frequency 

10 

kHz 

25 

kHz 

50 

kHz 

Re  '  Z,  ) 

Im'Z,  ) 

Re- Z,  ) 

Im'Z,  ) 

Re' Z,  ) 

Im '  Zi 

depth 

2  m 

84351 

40792 

104349 

23921 

109883 

15107 

4 

89755 

30084 

101470 

1  9481 

107955 

1  5003 

6 

87953 

27200 

100796 

20026 

108073 

1  51  09 

8 

86902 

27245 

100904 

201  21 

1 08070 

1  5097 

10 

86762 

27540 

10091  6 

201  00 

108070 

1  5097 

12 

86823 

27620 

100913 

201  00 

1 08070 

15097 

14 

86854 

27615 

10091  3 

201  00 

1 08070 

15097 

1  6 

86857 

27606 

100913 

20100 

108070 

1  5097 

1  8 

86855 

27604 

10091  3 

20100 

108070 

1  5097 

20 

86854 

27604 

10091  3 

20100 

1 08070 

1  5097 

22 

86854 

27604 

10091  3 

201  00 

108070 

1  5097 

Table  3*a  High-frequency  data.  Data  with  a  0.5  % 

random  error  are  subjected  to  nonlinear  least-squares 
inversion.  Compared  are  data  without  error  and  theore¬ 
tically  fitted  values. 

/r  =  h=  z=  10m,  (7,  =  2.0  ,  CT*  =0.2  mho/m,  d  =  12  m)  . 


data  with  0.5  $ 
error 


data  without 
error 


the  theoretically 
fitted 


25  kHz 

50  kHz 

25  kHz 

50  kHz 

25  kHz 

50  kHz 

Re  1 Z,  ) 

100999 

108263 

100913 

1 08070 

100895 

108051 

Im(Z,  ) 

20070 

15118 

20100 

1  5097 

20096 

15095 

Re( Z* ) 

-17419 

-17739 

-1 7428 

-17718 

-17426 

-17717 

Im(Zf ) 

-1  640 

-906 

-1642 

-904 

-1643 

-905 

Table  3 • b  Nonlinear  least-squares  inversion  results 
'  h  and  01  )  . 


Initial  value  Final 


h 

O', 


8.0  m 

1 . 6  mho/ m 


10.0012  m 

1 • 9997  mho/m 
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SELECTION  OF  OPERATING  FREQUENCIES 

Once  the  data  are  corrected  properly  for  the  misposi- 
tion  of  the  bird  through  the  methods  presented  in  the  previ¬ 
ous  sections,  the  next  problem  is  to  determine  which  data 
(in  frequency)  contain  pertinent  information  toward  a  suc¬ 
cessful  bathymetric  interpretation  for  a  given  range  of 
water  depth.  In  this  section,  first  we  shall  concentrate  on 
how  to  determine  an  optimum  frequency  range  for  a  given 
range  of  water  depth.  Next  we  will  see  an  example  of  inter¬ 
pretation  through  the  use  of  synthetic  data. 

Optimum  Frequency  Range 

One  of  the  first  task  for  an  airborne  EM  bathymetry 
survey  is  to  determine  the  optimum  frequency  range  for  a 
given  range  of  water  depth.  More  importantly,  what  is  the 
lowest  frequency  without  which  the  depth  determination 
becomes  unreliable?  We  shall  start  answering  this  question 
by  considering  the  following  example. 

Let  us  consider  inversions  of  eight  different  water 

depths  d  =  10,  20,  -  ,  80  m,  for  conductivities  of  sea 

water  =  2  mho/m  and  bottom  sediments  <Tt  =  0.2  mho/m. 

Data  are  generated  theoretically  for  six  frequencies  between 
50  Hz  and  15-8  kHz,  equally  spaced  on  a  logarithmic  scale 
(0.05,  0.16,  0.5,  1.6,  5  and  15.8  kHz),  at  a  10  m  bird  alti¬ 
tude  with  a  10  m  coil  separation.  For  a  theoretically  gen¬ 
erated  data  set  without  error,  we  find  that  the  nonlinear 
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inversion  scheme  yields  exact  answers.  For  a  data  set  with 
errors  (real  data  cannot  be  error-free),  the  scheme  can 
always  be  made  to  converge  (Aki  and  Richards,  1980). 

Table  4 • a  shows  the  frequencies,  skin  depth,  and  a  set 
of  contaminated  data  (for  d  =  10  m) .  Table  4-b  illustrates 
the  inverted  results  using  amplitude  of  coplanar  system. 

The  results  indicate  that,  i)  water  conductivity  is  deter¬ 
mined  well  within  0.25  %  regardless  of  depth  of  water,  ii) 
sediment  conductivity  is  not  determined  accurately.  More¬ 
over  after  50  m,  it  is  not  resolved  at  all  (in  the  inversion 
process  the  sediment  conductivity  is  allowed  to  vary  between 
0.4  -  0.001  mho/m),  iii)  uncertainty  in  the  depth  determina¬ 
tion  also  increases  along  with  that  in  the  sediment  conduc¬ 
tivity  for  depths  greater  than  50  m.  Error  bounds  in  the 
computed  depth  are  shown  in  parenthesis  in  Table  4*b.  The 
lower  bound  is  found  by  expanding  the  water  depth  until  the 
rms  difference  between  the  data  and  the  predicted  response 
becomes  more  than  0.25  %  assuming  an  insulating  ocean  bot¬ 
tom.  Similarly,  the  upper  bound  is  found  by  assuming  0.4 
mho/m  as  the  upper  limit  of  the  sediment  conductivity.  The 
width  of  the  bounds  starts  to  increase  rapidly  for  water 
depth  greater  than  50  m. 

As  expected,  the  above  results  confirm  the  fact  that 
determination  of  water  conductivity  is  not  problematic  for  a 
given  frequency  range  regardless  of  water  depth,  however, 
determination  of  water  depth  and  sediment  condcutivity  is 
highly  dependent  upon  the  lowest  frequency  used.  The  prob- 
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lem  of  determining  the  maximum  lowest  frequency,  for  the 
practical  purpose,  still  remains. 

In  the  linear  inverse  theory,  determination  of  the 
unknown  parameters  coincides  with  solving  a  linear  matrix 
equation: 

y  =  A  x  '24) 

where  A  is  an  nxm  coefficient  matrix  (n  is  number  of  data,  m 
is  number  of  unknown  parameters),  x  is  a  solution  vector  and 
y  is  a  vector  which  is  subjected  to  be  fitted  with  the  theo¬ 
retical  model.  In  the  generalized  inverse  theory,  the  theo¬ 
retical  curve  y  is  given  by 

y  =  S  y  '25) 

where  y  and  y,  respectively,  denote  the  predicted  and  the 
observed  and  S  is  the  information  density  matrix  'Jackson, 
1972;  Glen  et  al . ,  1973)*  We  note  that  y  is  a  weighted  sum 
of  y  through  the  rows  of  information  density  matrix.  Each 
element  in  a  row  of  S  represents  the  relative  contribution 
of  the  observed  to  the  predicted.  The  diagonality  of  S 
matrixs  indicates  whether  or  not  the  corresponding  datum 
significantly  affects  the  solution. 


Table  4 • a  Theoretically  generated  data  (d  =  10  m)  at 

six  frequencies  (0.5  -  15*8  kHz). 

(z  =  h  =  r  =  10m,  0)  =  2 ,  O'*.  =  0.2  mho/m). 


Frequency 
( kHz ) 


Skin  depth 

/  m) 


50 

50.3 

542 

1  58 

28.3 

3379 

500 

15-9 

1  561  0 

1  ,580 

8.95 

43502 

5,000 

5-03 

72682 

1 5 , 800 

2.83 

94424 

Data  with  0.5  %  error 
in  ppm  ( d  =  10  m) 

Real(Zi)  Imag.(Z|) 


4208 

11675 

25543 

35465 

31914 

23830 


Table  4-b 


(nonlinear  least-squares)  inverted  results 


true 

depth 

Cm) 


(T{  =  2  mho/ m 

cr; 

2.005 
2.004 
1  .999 
2.001 

Li-9,92 _ 

2.000 

1.993 

1.995 


01  =  0 . 2  mho/m 
Inverted  results 

01  d 

0.228  9-79 

0.179  20.10 

0.331  29.65 

0.213  38.20 

0. 1  96 _ 50.81 

0.027  61.70 

0.399  48.26 

4.000  48.88 


' bounds' 


( mho/m) 


The  information  density  matirices  for  each  depth  in  the 
example  are  shown  in  Figure  8.  The  bottom  rows  of  the 
matrix  recede  from  a  diagonal  shape  as  the  depth  increases. 
This  implies  that  the  high-frequency  data  becomes  less 
important  as  the  depth  increases.  On  the  other  hand,  the 
contibution  of  lower  frequencies  increases  as  the  depth 
increases . 

Examination  of  the  information  density  matrix  reveals 
that  i)  the  higher  frequency  data  decrease  their  contribu¬ 
tion  as  depth  increases.  There  is  a  certain  depth  where  the 
information  density  of  some  high-frequency  data  cease  to  be 
diagonal  (shown  by  dotted  line  in  the  information  density 
matrix) ,  and  that  part  of  the  data  do  not  contribute  in  con¬ 
struction  of  a  solution,  and  ii)  the  low  frequency  data 
increase  their  contribution  as  the  depth  increases;  i.e., 
there  is  a  certain  depth  where  the  diagonal  element  of  a  low 
frequency  becomes  unity  ('saturated').  Below  that  depth, 
the  contribution  of  the  data  is  saturated  and  for  a  reliable 
result  of  inversion,  lower  frequency  data  are  needed.  This 
may  be  interpreted  as  follows;  The  coefficient  matrix  A  in 
equation  (24)  consists  of  the  derivatives  with  respect  to 
model  parameters,  and  the  information  density  matrix  is 
defined  by  the  data  space  eigenvectors  of  A  (Lanczos,  1961; 
Aki  and  Richards,  1980).  Derivatives  above  a  certain  high 
frequency  become  much  smaller  than  those  of  lower  frequen¬ 
cies.  Those  vanishing  derivatives,  therefore,  do  not  con¬ 
tribute  in  solving  the  equation  '24).  Instead,  they  result 


in  the  nondiagonality  of  the  information  density  matrices. 
Therefore,  those  high  frequency  data  can  be  discarded  safely 
(e.g.,  dotted  line  in  the  information  density  matrix  is  the 
starting  point  of  nondiagonality) .  Meanwhile,  the  weight  of 
derivatives  at  lower  frequency  increases  as  the  depth  of 
water  increases  because  then  high-frequency  components  in  a 
column  of  'A'  become  similar.  At  a  certain  point  the  weight 
of  the  low-frequency  part  is  saturated.  At  that  point  it  is 
necessary  to  have  another  lower  frequency  to  achieve  satis¬ 
factory  resolution.  In  this  example,  the  depth  of  satura¬ 
tion  of  the  lowest  frequency  (dashed  line  in  S  plots) 
roughly  coincides  with  the  skin  depth  of  that  frequency 
(frequency  and  skin  depth  are  shown  at  the  bottom  of  each 
figure),  and  with  the  underlined  depth  of  inverted  result. 
The  deviation  from  diagonality  of  S  occurs  when  the  depth  of 
water  exceeds  3-4  times  of  the  skin  depth.  Similar  analy¬ 
ses  with  different  sets  of  frequency  contents  (Son,  1985)  do 
not  violate  the  above  idea  that,  if  the  depth  of  investiga¬ 
tion  is  deeper  than  the  skin  depth  of  the  lowest  frequency 
of  the  data,  we  need  lower  frequency  data  to  achieve  proper 


resolution. 
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Interpretation  Example 

Figure  9  illustrates  an  interpreted  bathymetry  profile 
based  on  a  set  of  simulated  AEMB  data  using  six  frequencies 
between  50  Hz  and  15*8  kHz,  equally  spaced  on  a  logarithmic 
scale  (0.05,  0.16,  0.5,  1.6,  5.  15-8  kHz).  Data  are  gener¬ 
ated  from  a  two-layer  model  in  which  the  first  layer  is  sea 
water  having  a  conductivity  of  2  mho/m  underlain  by  a  uni¬ 
form  bottom  sediment  having  a  conductivity  of  0.2  mho/m. 
Simulated  AEMB  data  are  obtained  for  forty  different  water 
depths  between  0  to  80  m  at  a  2  m  interval  at  a  10  a  alti¬ 
tude  with  a  10  m  coil  separation.  A  set  of  random  numbers, 
whose  rms  amplitude  is  normalized  to  that  of  an  infinitely 
long  sequence,  is  added  to  contaminate  the  theoretical  data 
by  0.5  £.  Dotted  lines  in  Figure  9  show  the  inverted 
depths,  while  the  solid  lines  represent  true  depths.  In  all 
40  inversions,  the  water  conductivity  is  determined  well 
within  0.25  '%  (1-9975  -  2.0046  mho/m). 

The  uncertainties  in  the  water  depth  and  sediment  con¬ 
ductivity  increase  after  about  40  m.  The  solid  vertical 
error  bars  represent  the  uncertainty  in  depth.  The  lower 
bound  is  found  by  expanding  the  water  depth  until  the  rms 
difference  between  the  data  and  the  predicted  response 
becomes  more  than  0.25  assuming  an  insulating  ocean  bot¬ 
tom.  Similarly,  the  upper  bound  is  found  by  assuming  0.4 
mho/m  as  the  conductivity  of  the  bottom  sediment.  The  Edge 
model  error  bounds  start  increasing  rapidly  after  50  m. 


reaching  infinity-  after  70  m. 

The  information  distribution  (diagonal  elements  of 
information  density  matrix)  is  illustrated  in  the  inset .  It 
shows  that  the  contribution  from  low-frequency  data 
increases  as  depth  increases,  while  the  high-frequency  con¬ 
tribution  increases  as  the  water  depth  becomes  shallow.  The 
optimal  frequency  range  (where  the  amplitude  of  diagonal 
elements  of  the  information  density  matrix  become  bigger 
than  1/2)  as  a  function  of  water  depth  is  shown  by  a  shaded 
area.  The  frequency  bands  within  this  area  contribute  most 
to  the  solution. 

CPU  times  for  the  first  ( d  =  2  m)  and  the  last  ( d  =  80 
m)  inversion  are  0.86  sec  and  3*27  sec,  respectively,  on  an 
IBM  3081  mainframe  computer.  All  least-squares  inversion 
methods  require  a  set  of  initial  values  which  is  subseq¬ 
uently  improved  through  iterations  (90  #  of  the  model  was 
used  as  an  initial  guess  in  this  inversion) .  The  required 
computer  time  is  determined  mainly  by  the  number  of  itera¬ 
tions  needed  to  achieve  a  prescribed  resolution.  In  order 
to  reduce  the  number  of  iterations,  we  have  devised  a 
"tracking"  scheme  in  which  the  initial  solution  of  the  pres¬ 
ent  location  data  will  be  prescribed  as  the  final  solution 
for  the  last  location  data  point.  This  scheme  is  ;ustified 
as  long  as  the  bathymetry  changes  gradually  along  the  flight 
path.  If  successful,  the  scheme  can  promote  significant 
saving  in  computer  time  for  AEMB  data  processing. 


Figure  9  Frequency-domain  EM  data  inversion  simulation, 

showing  the  bathymetry  resolution  up  to  80  m  of  water  depth 
Simulated  data  are  contaminated  by  0-5  %  random  measure¬ 
ment  errors  for  all  cases.  The  sea  water  conductivity  is 
resolved  within  0.25  #  in  all  cases.  Resolution  of  the 
sediment  conductivity  is  shown  at  the  bottom.  The  freq¬ 
uency  bands  which  contribute  most  to  the  inversion  are 
shown  as  a  shaded  area  in  the  inset. 


Bothymetry  FREQUENCY- DOMAIN  EM  DATA  INVERSION  SIMULATION 

_  WATER  SURFACE 


BOTTOM  CONDUCTIVITY  AS  A  FUNCTION 


EXPONENTIALLY  DECAYING  CONDUCTIVITY  MODEL  UNDERLAIN 


BY  SEA  WATER 

Because  of  high  conductivity  of  sea  water,  the  conduc¬ 
tivity  of  sediment  is  expected  to  be  proportional  to  the 
porosity  of  bottom  sediment  rather  than  the  conductivity  of 
sediment  material.  Then  the  conductivity  in  the  sediment 
layer  will  be  decaying  with  depth,  rather  continuously  than 
abruptly  like  in  the  layered  earth  model. 

The  continuously  changing  conductivity  models  appear  in 
a  few  literatures.  Mallick  and  Roy  (1971)  has  derived  a 
formula  for  the  response  of  vertical  magnetic  dipole  over 
transitional  earth  where  the  linear  transition  zone  exists 
between  the  overburden  and  the  lower  half-space.  Puller  and 
Wait  (1972)  has  derived  an  expression  for  the  high-frequency 
electromagnetic  coupling  between  small  coplanar  loops  over 
an  inhomogeneous  ground,  where  the  conductivity  of  tran¬ 
sition  zone  under  resistive  overburden  increases  exponen¬ 
tially  toward  the  larger  values  of  O^deep  within  the  earth. 
In  this  paper,  we  present  a  theoretical  electromagnetic 
dipole  response  over  exponentially  decaying  conductive 
material  underlain  by  sea  water,  with  the  aim  of  investigat¬ 
ing  the  conductivity  of  bottom  sediment  of  sea  water. 

The  conductivity  model  is  shown  in  Figure  10. a.  The  (Jj 
and  d  represent  the  conductivity  and  thickness  of  the  sea 
water,  respectively.  The  conductivity  O^of  the  bottom  sedi¬ 
ment  is  given  by 


where  q  represents  the  decaying  factor,  and  z  represents  the 
depth.  In  this  model,  since  z  is  less  than  -d,  z+d  becomes 
negative.  Therefore,  if  q  is  big,  (^decays  fast,  if  q  is 
small,  (Ji  decays  slowly.  The  conductivity  of  bottom  layer  is 
equal  to  that  of  sea  water  at  the  interface.  At  1/q  m  below 
the  interface,  it  becomes  1/e  of  the  conductivity  of  sea 
water.  Eventually,  it  becomes  zero  deep  in  the  earth. 

The  kernel  function  R'A)  in  expression  MO)  can  be 
replaced  by  (27)  for  the  exponentially  decaying  conductivity 
model . 


R'A) 


•n.-x, 


x.  -n, 


+  -0/(1 


X.  -  X, 


- JL  *.  ) 

n.  -*  -vi, 


.'2  7) 


where  =  - — • — --  eZn,<i  (see  Appendix), 

p-n,  -»  p' 


When  d  is  set  to  infinity  in  equation  (27),  becomes  zero 

and  R(A)  =  - - ,  which  is  the  case  of  homogeneous  half- 

space.  As  q  approaches  infinity,  the  conductivity  of  sedi¬ 
ment  decays  rapidly  and  the  model  becomes  the  case  of  two- 

layered  earth  where  the  sea  water  layer  is  overlain  by 

"w  “  A 

insulating  bottom  layer.  In  this  case,  <*.  approaches  — - - e 

I'll  +  A 

and  equation  '27)  becomes  identical  with  equation  (10). 

Figure  10. b  shows  responses  of  the  exponentially  decay¬ 


ing  conductivity  model  with  varying  q  values.  In  Figure 
10. c,  compared  are  the  responses  of  two  layer  model  with 
varying  bottom  layer  conductivities.  The  responses  are  cal- 


culated  at  42.67  m  (140  feet)  altitude  with  8  m  coil  separa¬ 
tion.  Inphase  component  increases  as  (^increases  or  q 
decreases.  The  rate  of  increase  in  exponentially  decaying 
conductivity  model  is  smaller  than  that  in  two  layer  model, 
especially  at  lower  frequencies,  indicating  the  conductivity 
of  the  former  model  is  smaller  than  that  of  the  latter  at 
below  a  certain  depth.  The  peak  values  of  quadrature  compo¬ 
nents  appear  at  lower  frequencies  and  becomes  smaller  as  CT*. 
increases.  The  same  phenomenon  happens  as  q  decreases, 
indicating  that  eddy  current  is  spread  in  a  wide  range.  The 
decay  rate  of  peak  amplitude  in  this  case  is  smaller  than  in 
two  layer  model. 

Figure  11  compares  the  responses  of  four  different  mod¬ 
els.  Model  1  in  Figure  11. a  is  a  two  layer  model  where  the 
conductivity  and  the  depth  of  sea  water  are  4  mho/m  and  4  m, 
respectively,  and  the  conductivity  of  bottom  layer  is  0.01 
mho/m.  Model  2  is  that  of  half-space  with  conductivity  of  4 
mho/m.  Model  3  is  an  exponentially  decaying  conductivity 
model  with  q  =  1,  underlain  by  4  m  thick  sea  water.  In 
model  4,  the  exponentially  decaying  part  of  model  3  is 
devided  into  fifty  horizontal  layers  each  having  constant 
thickness  of  10  cm.  The  responses  are  plotted  in  Figure 
1 1 . b  (the  response  of  model  4  is  represented  by  '+')•  The 
responses  of  model  3  and  model  4  are  in  agreement. 


Figure  10  The  exponentially  decaying  conductivity  model* 

a.  The  exontially  decaying  conductivity  model  underlain 
by  4  m  thick  sea  water.  The  01  and  Qi  represent 
conductivity  of  sea  water  and  sediment,  respectively, 
and  d  represents  the  depth  of  water  layer. 

b.  Response  curves  for  the  exponentially  decaying  conduc¬ 
tivity  model  with  different  q  values  MOO,  2,  1,  1/2, 
1/4)  . 

c.  Response  curves  for  two  layer  model  with  different  0 

( 2 ,  1,  0.5,  0.2,  0.001).  The  responses  are  calculated 
at  140  feet  (42.67  m)  with  8  m  dipole  separation. 


Figure  11  Comparison  of  responses  of  four  different  models. 

a.  Geometry  of  four  different  models.  Model  1  is  a  two 
layer  model  where  the  conductivity  and  the  depth  of 
sea  water  are  4  mho/m  and  4  m,  respectively,  and  the 
conductivity  of  bottom  layer  is  0.01  mho/m.  Model  2 
is  that  of  half-space  with  conductivity  of  4  mho/m. 
Model  3  is  an  exponentially  decaying  condcutivity 
model  with  q  =  1,  underlain  by  4  m  thick  sea  water. 

In  model  4,  the  exponentially  decaying  part  of  model 
3  is  devided  into  fifty  horizontal  layers  each  having 
constant  thickness  of  10  cm. 

b.  The  responses  of  four  different  models.  The  numbers 
on  the  curve  coincide  with  those  of  models.  The 
response  of  model  4  is  plotted  by  '+’. 


CONCLUDING  REMARKS 


In  the  Airborne  Electromagnetic  Bathymetry  problem,  the 
ma^'or  challenge  is  the  skin  effect  of  highly  conductive  sea 
water.  Since  most  of  the  response  is  due  to  uppermost  sea 
water,  errors  in  sensor  altitude  and  inclination  can  cause 
as  much  effect  as  bottom  geometry.  To  increase  water  depth 
sensitivity,  we  are  forced  to  lower  the  frequency  band, 
which  also  suffers  from  severe  practical  limitations. 

Despite  all  these  difficulties,  the  AEMB  concept  carries  a 
few  advantageous  aspects  including  a)  a  simple  geometry 
which  may  be  approximated  as  one  layer  over  a  homogeneous 
half-space,  b)  no  surface  topography,  c)  a  more  or  less  uni¬ 
form  sea  water  conductivity  within  a  zone  below  the  instru¬ 
ments  (due  to  rapid  mixing  of  shelf  water  caused  by  wave 
action) ,  and  d)  possible  application  of  bathymetric  tracking 
during  the  inversion  procedure  owing  to  gentle  and  gradual 
changes  in  bathymetry.  These  simple  geometrical  aspects 
provide  us  with  many  practical  advantages  needed  in  survey 
design  and  data  interpretation. 

The  results  in  this  study  include  i)  that  uncertainties 
in  sensor  altitude  can  be  resolved  through  an  inversion  pro¬ 
cess,  using  high-frequency  data  where  the  effect  of  bottom 
sediments  is  excluded,  and  assuming  two  unknowns  (altitude 
and  sea  water  conductivity),  ii)  inclination  of  the  coil 
axis  with  respect  to  the  sea  surface  can  be  resolved  by 
installing  an  inclinometer  on  the  sensor  platform  and  cor- 


recting  for  this  effect  during  the  inversion  process,  iii) 
the  optimum  frequency  range  for  operational  and  interpreta- 
tional  purposes  can  he  determined  before  actual  measure¬ 
ments,  via  information  density  matrix  analysis.  For  the 
data  measured  at  10  m  altitude  with  10  m  separation,  it 
appears  that,  for  a  reliable  result,  at  least  a  low  fre¬ 
quency  whose  skin  depth  is  close  to  the  depth  of  investiga¬ 
tion  need  be  included,  and  that  high  frequencies  whose  skin 
depths  are  less  than  1/4  -  1/3  of  the  depth  of  investigation 
are  not  important  in  determining  the  bottom  geometry. 
Interpretation  and  survey  design  may  vary  depending  with 
each  situation  (percentage  of  error  included  in  the  data, 
height  and  separation  of  the  bird,  etc.).  However,  the 
methods  presented  in  this  study  could  be  applied  in  any 
case.  Vertical  heterogeneity  of  sea  water  or  bottom  sedi¬ 
ment  conductivity,  and  deviation  from  a  horizontally  layered 
model  were  not  considered  in  the  analysis. 

In  order  to  account  for  the  effect  of  conductivity  var¬ 
iation  in  sediment,  a  theoretical  electromagnetic  response 
over  exponentially  decaying  conductivity  material  underlain 
by  sea  water  has  been  derived.  Comparison  of  the  response 
with  that  of  a  similar  multi-layered  model  shows  agreement. 
Since  the  decaying  pattern  of  sediment  conductivity  is 
related  with  the  compactness  of  the  bottom  sediment,  it  may 
provide  valuable  information  toward  ocean  engineering. 
Application  of  the  theory  presented  in  this  paper,  to  real 
data,  will  be  published  in  next  paper. 


APPENDIX 


Derivation  of  a  kernel  function  for  exponentially- 
decaying  conductivity  model 


In  this  appendix,  we  will  derive  a  response  function 
(mutual  coupling  ratio)  of  the  horizontal  coplanar  system 
over  a  exponentially  decaying  conductivity  model  underlain 
by  highly  conductive  sea  water  (Figure  10. a).  The  electri¬ 
cal  permittivity  («)  and  the  magnetic  permeability  (m)  of 
the  earth  will  be  assumed  the  same  as  those  of  the  air.  The 
displacement  current  in  the  earth  is  neglected. 

In  time-harmonic  electromagnetic  field,  the  magnetic 
vector  potential  A  satisfies  the  Helmholtz  equation  (Patra 
and  Mallick,  1980): 


z  -*  ,  -> 

V  A  +  Y1  A  =  0 


( A-1  ) 


where  =  iw)u<5\  i  =  V -1  ,  w  angular  frequency  and 
G  is  conductivity  of  the  material.  For  the  horizontal 
dipole  over  vertical  stratification,  the  vector  wave  equa¬ 
tion  (A-1)  reduces  to  a  scalar  wave  equation  in  cylinderical 
coordinate  (r,  9,  f  ) : 


+  1  $-1  , 

"i"rl  r  Tr  “?r 


+  i*A§  = 


o 


'A-2) 


where  Ax  is  z  component  of  the  vector  potential  A. 
Separating  the  variables:  A(=Az)  =  R(r)  Zrz),  and  substi¬ 
tuting  this  expression  into  equation  (A-2),  we  obtain  two 
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ordinary  differential  equations: 


+  -  -~  +  A2,  R  =  0 

(A-3) 

"oTT1- 

It  dr 

d2Z 

~6  il 

-  {  *l(z)  +  A2  |  Z  = 

0 

(A-4) 

where  A  is  a  separation  constant.  A  non-divergent  solution 
of  equation  (A-3)  is  J0 (Ar) :  Bessel's  function  of  order 
zero.  The  solution  for  equation  (A-4)  depends  upon  Tslf  z)  . 
For  a  layered  earth  model,  a  general  solution  of  the  vector 
potential  above  the  surface  of  the  earth  is  written  as  a 
Hankel  transform  integral: 

A  =  C  {  i  +  (  R'a)  e'A<k~*;  Ja  (Ar)  d*  J  (A-5) 

where  C  is  a  constant,  r  is  the  distance  between  two  loops, 
h  and  z  are,  respectively,  the  heights  of  transmitting  and 
receiving  loops,  and  R(A)  is  a  kernel  function  including 
geological  parameters.  The  kernel  function  for  a  two  lay¬ 
ered  earth  is  given  by  equation  (10). 

When  we  consider  the  exponentially  decaying  conductiv¬ 
ity  model  underlain  by  highly  conductive  sea  water  (Figure 
10. a),  the  radiation  constants  in  each  region  can  be  written 
as , 


(sea  water)  : 

(sediment)  : 

n* 

iw/iu; 
iWyU  C(  z) 

»,*  e*'1^ 


Then,  the  equation  (A-4)  in  the  sediment  layer  becomes, 


f*  e*c*+Jj  +  A1  ) 


Making  a  substitution  of  the  form 


x(z) 


we  can  derive 


-  d2Z 
and  ---- 
a?* 


d*Z  f  <*•  JZ.  X 
(q  x)  +  T-  q  x 


With  these  expressions,  equation  (A-6)  becomes 

x*  7?-  +  q1  x  #-  -  (x  +  a1)  %  =  o 

In  this  equation  if  q  is  not  equal  to  zero  (If  q  =  0,  it 
becomes  the  case  of  two  layer  model.)*  we  obtain  a  linear 
homogeneous  second-order  differential  equation  with  a  regu 
lar  singular  point  (Arfken,  1971)  at  x  =  0: 

+  \  JZ  x  +  A*  z  _  o  'A_ 

Tx1  a  4*  {***■ 

A  solution  of  this  equation  can  be  secured  by  Probenius' 
method  assuming  the  solution  of  the  form, 


Z(x) 


a  x 
i 


(a.  j  0) 


where  k  is  a  constant.  Substituting  the  expression  (A-8) 
into  equation  (A-7),  we  derive  a  non-diverging  solution 


Z(x) 


-  x(*  "  ^ 

C  — 

a*0  0*)/ 


where  <  x  =  r,*e*  '<  0,  q.  ^  0  and  (  +  2\)  t  =  1  ,  for  r1  =  0. 


Let  us  write 

.iL-'-ll-l-}------ 

(4%) !  (i%  »a;  / 

Then,  Z(z)  =  a0  P(z) 

The  corresponding  solutions  in  different  layers  are: 


*  =  • 


Air  :  A  =  C  i  -p  +  |  A(\)e'1'*Z  J„  (Ar)  dA  } 

Sea-Water  :  A*  =  (  { B ,  ( A)  e""'Z  +  B2(a)  e*,Z  } 

Lower  Half-Space  :  A3  =  j  D'A)  P(z)  J0 'at)  dA 

where  nt  =  /  a1  +  v*  ,  i  =  0,  1,  2.  The  A(a),  B,  (A)  ,  Bz(a) 
and  D(a)  are  coefficents  which  will  he  determined  from 
boundary  conditions. 

Applying  the  boundary  conditions  (the  tangential  compo¬ 
nents  of  elctric  and  magnetic  field  vectors  must  be  continu¬ 
ous  accross  the  surface)  at  z  =  0  and  -d,  we  can  derive 


A(  A ) 

t  Tio’-n. 

V  -n.  ♦  -n, 

Pn,  -  p' 

P  "»»i  ■*  P' 

ao 

p 

p(-d)  =  r 

and  P  ’ 

=  5SL\  = 

+  -0/(1  +  « 


n, 


-  i-n.  i 


V. 


£  Lt'j-V. 

*  =  '  (*Z>!  <*%  +  **)! 


(A-10 


This  is  the  kernel  function  for  exponentionally  decaying 
conductivity  model  underlain  by  sea  water. 


Convergence  Consideration 


Before  we  consider  numerical  convergence,  let  us  con¬ 
sider  two  limiting  case  of  the  kernel  function.  When  d  is 
set  to  infinity  in  equation  (A-10),  the  model  becomes  a 
half-space  having  sea  water  conductivity.  In  this  case  <*- 
becomes  zero  and  R(A)  =  which  is  identical  with  the 

'H.-Wl, 

kernel  function  of  a  half-space.  This  is  a  confirmation  of 
the  result  in  (A-10).  When  q  becomes  infinity,  the  model 


approaches  a  two-layered  one  having  zero  conductivity  in  the 
second  layer.  In  this  case,  P  and  P'  in  (A-10)  become,  lim 

%-9°° 

P  =  1  and  lim  P'  =  A  .  Then, 


lim 


'vi-  A  -i  n,  d 

- e 

y\.  A 


This  being  substituted  into  (A-10),  the  kernel  function 
becomes  identical  with  expression  (10)  for  a  two  layer  model 
having  zero  second  layer  conductivity.  This  is  another  con¬ 
firmation  of  the  result. 

Convergence  of  P  and  P'  in  (A-10)  is  dependent  upon  the 
value  of  q.  If  q  is  big  they  converge  fast.  If  q  is  small 
they  converge  slowly.  The  decaying  factor  q  determines  the 
decay  rate  of  sediment  conductivity.  For  example,  if  q  =  1 
(1/m),  (Jz  becomes  1  / e  of  0)  at  1  m  below  the  water-sediment 
interface.  If  q  =  0.1,  0^  becomes  1/e  of  0(  at  10  m  below 
the  interface.  At  q  =  oo  which  is  the  case  of  two  layer 
model  having  infinitively  resistive  second  layer,  the  first 
terms  in  the  series  are  dominant.  As  q  decreases  (01  decays 


slowly) *  the  higer  terms  are  added  and  convergence  becomes 
slower . 

Referring  back  to  the  derivation  of  equation  (A-7),  the 
result  in  (A-10)  is  valid  only  for  q  =  0.  Moreover  as  q 
approaches  zero,  numerical  convergence  may  confront  serious 
problem  (limitation  of  the  result) .  However,  the  case  of 
extremely  small  q  (less  than  0.01)  is  rare  in  reality.  If 
we  are  considering  q  bigger  than  0.1  like  in  the  bathymetry 
problem,  the  convergence  is  not  problematic. 
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